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A VORTEX MODEL FOR TRANSPORT IN THE 
POLAR STRATOSPHERE 

I. J. Eberstein, F. Y. Yap and V. Viers 
ABSTRACT 

A semi-empirical model based on a Gaussian vorticity distribution has been developed 
for determining eddy diffusivity and wind transport distributions in the polar stratosphere. 
The model uses as input data pressure surface heights measured at periods of the year when 
the stratospheric polar vortex, exhibits nearly circular patterns around the pole. The com- 
ponents of the polar wind velocities that result from a Prandtl eddy viscosity distribution 
are found to be in general agreement with those obtained by other investigators. 
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A VORTEX MODEL FOR TRANSPORT IN THE 
POLAR STRATOSPHERE 

A. INTRODUCTION 

A vortex based model has been developed to calculate the eddy diffusivity distribu- 
tion and the components of the wind velocities in the polar stratosphere and mesosphere 
in the altitude range 1 5-60 km. The model is semi-empirical in its application, and re- 
quires input data of measured pressure surface heights between lOOmb and 0.4mb during 
the late fall-early winter months when the pressure contours exhibit well developed and 
relatively undistorted circular patterns centered at or near the north pole. In addition to 
determining the wind transport, the model is also useful in investigating the relationship 
between eddy viscosity and mean flow, and the effect of these parameters on polar 
stratospheric ozone. 

This paper describes the mathematical development of the polar vortex model, and 
the numerical computer techniques that are employed in performing the calculations. 
Results generated by the model are presented and discussed. 

B. THE EQUATIONS OF MOTION 

The geometric configuration and coordinate system that are used to develop the polar 
vortex wine 1 model are shown in Figures 1 and 2. The coordinate system is fixed to a 
point of observation O at latitude X on the surface of,the Earth which is rotating with 
angular velocity <3. In this rotating coordinate system the velocity vector ^ of a particle 
may be written in cylindrical coordinates (r, 0, z) as 
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Figure 1. Geometric Configuration of a Rotating Coordinate 
System Fixed to the Earth’s Surface 
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Figure 2. Cylindrical Coordinate System Used in the Development 
of the Polar Vortex Wind Model 
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( 1 ) 


Vr ^ , a /s 

V = vr + up + wz , 


where 

'r', $> and 'z are unit vectors in the r, 4> and z directions, respectively, 
u = the zonal component of velocity, 
v = the meridional component of velocity, 
w = the vertical component of velocity. 

We note that at large latitudes X when O is close to the north pole, the curved distance (ft 
from the pole to O closely coincides with the r axis of the rotating cylindrical coordinate 
system. 


The general equation of motion for an element of unit mass, with respect to the 
rotating coordinate system is (Craig, 1960): 

•> d V l » 

a = _ = _ _ VP _ 2c$ x V + g + f , (2) 

dt p 


where 

a is the -acceleration of the unit mass of air, 

— VP is the force due to the pressure gradient, 

P 

- 2(co x \/) is the Coriolis force due to the Earth’s rotation, 
g is the force due to gravity, which acts only in z direction, 
f is the frictional force, 
p is the air density. 
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We may express certain terms in the above equation in cylindrical coordinates as follows: 

a = (r - r^ 2 )'? + ( r0 + 2r 0)^> + zz 


-> 3P i 3P^ ap^ 

VP = — ? + - — $ + —z 

9r r 90 9z 


oj x V = - ucosinX? + (vtosinX + wcocosX)^ - ucocosXz' 
Therefore, Equation (2), in component form, becomes: 

1 3P 

r - r 0 2 + 2uccsinX + f 

p 9r 

1 dP 

r 0 + 2r0 = - — 2(vcosinX + wcocosX) + 1 'a 

pr 90 

1 9P 

z + 2ucocosX + g + f_ . 

p 9z 


(3a) 


(3b) 


(3c) 


C. THE ZONAL WIND SPEED 

The vortex model assumes a vorticity £ which has a Gaussian distribution with axial 
symmetry about the center of the coordinate system near the north pole: 


? = to e 


2 o‘ 


(4) 


where r is the radial distance from the center and a is the width of the distribution. 

^ 

From the definition of vorticity ? = V x V, Stokes’ Theorem may be applied around a 

o 

cylindrical contour of constant zonal speed: 




V x Vds = 


x 2 

/ * 0 e 2 ° 2 27rrdr = (j>V • dl 


= 2ttRu 


( 5 ) 


where u is the zonal component of the wind velocity and R is the radial distance from the 
center at the point where the zonal speed is to be calculated. 


Integrating Equation (5) we obtain 



We may define the vortex strength K and vortex radius A as 



(7) 


and obtain the zonal speed u as a function of K, A, and R: 

2K 


u = 


R 


1 - e 2A 




For the calculation of the zonal speed, we assume that the atmosphere has a mean 
constant zonal flow in circular patterns moving with angular speed fi relative to the 
center of the rotating coordinate system. This motion should not be confused with the 
angular rotation of the Earth. To a first approximation we may neglect the effect ; f 
friction on the zonal speed due to the weak radial gradients of u. Also, since r = constant, 
r = 0. The <p in Equation (3a) may also be written as 


d <j> _ ds d(j> _ u 
dt dt ds r 
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Hence tlie meridional equation of motion Equation (3a) becomes 


u 2 


1 3P 
p 3r 


+ 2ucosinX 


or, 


3P u 2 

— = p — ± 2ucesinX 
3r r 


( 9 ) 


The ± sign is introduced to differentiate between zonal flow toward the East (cyclonic) 
And flow toward the West (anticyclonic), respectively, in the Northern Hemisphere. 

( ( R \ / r \ 

Noting that sinX = cos I, — 1 — cos I — where R E is the Earth’s radius, we may 


Rr 


Rf 


integrate the above equation: 


r f R r » 2 / r \i 

I dP = / p — ± 2pucocos j — ) 
-'Pref ^RiefL f \ R e/_ 


dr 


( 10 ) 


where P re f is a chosen reference pressure and R re f is the distance from the pole where 
P re f is measured. R in the distance from the pole of the point at which the zonal speed 
u is to be calculated. 


The hydrostatic equation may be written: 


I dP = -g / pdz = - pg(z - z ref ) 
Jv r J? 


(H) 


"ref 


where we have assumed that p is equal to the average density over the pressure surface, 
Thus Equation (10) becomes 

1 f K f u 2 (R,K, A) 

l ref 


z z ref 


i r K r u 2 (r,i 

g L r 


2ucocos 


5)1 


dr . 


( 12 ) 
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where (R re f, z re f) defines the position of the reference point, chosen to correspond to the 
position of maximum latitude on the measured isobaric surface (see Fig. 3). (R, z) are 
the coordinates of the point at which the zonal speed u is to be determined. 

From Equation (8) we recall that u is a function of R, K and A. Equation (8) may 
be substituted into Equation (12) which is then solved to furnish the vortex strength K 
and vortex radius A as functions of the altitude z and the radial distance R from the 
origin. 

If z and R are known on an isobaric surface, then Equation (12) may be solved nu- 
merically for K and A. Using the following notations: 
n = index for a pressure surface, 
m = number of data points on the n^ pressure surface, 
i = index of data point on a pressure surface, 
we may rewrite Equations (8) and (12) for the n th contour as: 



From pressure contour maps published by the Environmental Science Sendees Administra- 
tion (ESSA, 1969, 1970), we obtain data points R re f (n , Zj >n , Rj n (i = 1 m) for the 

n^ 1 pressure contour. Typically, the pressure contours consist of isobars at 0.4, 2, 5, 10, 
30, 50, and 100 millibars. 
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O 

Figure 3. Illustration of Points on an Isobaric Surface 
Used in the Evaluation of Equation (12) 



Equation (13) is used to calculate U, n at each Rj n for the n ,h contour. For each 
contour there are m values of u, which are then used as input data to Equation (14). Note 
that K n and A„ are functions only of the particular contour, and do not depend on the 
distance from the pole. From the resulting in equations the best values of K n and A n can 
be found by the method of least squares for the n 1 *' contour. 


Since each contour may be considered as data taken at a certain average height, the 
K n ’s and A n ’s are therefore functions of height; curves may be fitted to them by least 
squares to provide K(z) and A(z) as continuous functions of height z within the altitude 
range of the data. From these cuives the zonal wind may be calculated as a function of 
latitude R and height z according to Equation (8): 


u(R,z) 


2K(z) 

R 


R 2 

2A(z)2 


1). Till: MERIDIONAL WIND SPEED 


(15) 


The effect of frictional forces on the meridional wind speed (which is roughly two 
orders of magnitude smaller than the zonal wind) is not negligible, as was the case for the 
calculation of the zonal wind. We rewrite the 0 equation of motion Equation (3b): 

1 91 ’ 

r 0 + 2r0 = - — 2(va>sinX + wcjcosX) + f0 (3b) 

pr 30 


, • v 3P 

Also, it we assume steady state motion and axial symmetry, r = 0, 0 = 0 and — = 0. 

30 

Hence the above equation becomes 


2(vtosinX + wcocosX) = L, 


(16) 
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Hie f rictional forces in the atmosphere are generally thought to be similar in behavior to 


forces due to eddy stresses arising from turbulent motion. It can be shown that the eddy 
stresses give rise to a force, the 0 component of which can be represented as (see Craig. 
I960): 


r 0 


I d 
p dy 



(17) 


where t is called the coefficient of eddy viscosity. 


The functional form of e is not completely known. Various forms have been 
suggested, among which is the following due to Prandtl (Hess. 1959): 

e = pi 2 1 — I s pe K (18) 

07 . 

where I is a characteristic mixing length, which is analogous to the molecular mean free 
path in a turbulence free situation, and c K is called the kinematic eddy viscosity. 


Since the vertical wind speed w is small compared to the meridional wind speed v. 
at points near the pole the w cocos X term in Equation (16) may be neglected, yielding: 


2wvsinX — t^ 


whence 



2pu>sin\ 


(19) 



Using this expression, we calculated the meridional wind v at any given location from the 
du 

zonal wind shear — . since the latter was readilv obtained from Equation (15) using nu- 
3z 

merical computer techniques. We note that the approximation in deriving Equation (19) 
improves at larger latitudes. 


E. THE VFRT1CAL WIND SPEED 

For calculating the vertical wind speed w, we make use of the equation of 
continuity 


3 p -» -t 
— + V • pV = 0 
3t 


(20) 


3 p 

Furthermore, if we consider the atmosphere to be in a steady state, — = 0. and the above 

3t 

equation reduces to: 


V • (pV) = 0 . 


( 21 ) 


In cylindrical coordinates (r, 0, z) the divergence of a vector A is given by: 


-► 1 3 I 3A, 

* A * - — (r A r ) + - + 


3A, 


r 3r 


r 30 3z 


Fquation ^2 1 ) may be rewritten as: 


tZ , -+ 1 3 1 3(pu) 3(pw) 

V • ( P V) = - - (prv) + - = 0 

r dr r 30 3z 


( 22 ) 


From axial symmetry — (pu) =0, and Equation (22) becomes 
30 


I 3 3 

- — + — (pwi = 0 

r 3r 3z 
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or. 


v dv dp dp dw 

p — + p— + v — + w — + p — = 0. 

r dr dr dz dz 


We assume that the equation for hydrostatic equilibrium is valid: 


dP = - gpdz . 


Using the ideal gas law and assuming isothermal layers (i.e., the temperature is approxi- 
mately constant over the thin layer between /.j and z 2 ). we may integrate the above 
equation to yield: 


= (z 2 - z, ) 

RT ‘ 


where P, and P 2 are pressures measured at the heights /., and z 2 . respectively. We may 


therefore write: 


]_ dp I 1 UP 2 /P 1 ) 

p dz z 2 - z\ 


which is substituted into Equation (23) to give 


dw v dv v dp ln(P 2 /Pi) 

dz r dr p dr z 2 - z, 


Since the density change with respect to latitude is small, the — term may be neglected. 

dr 

0 . 

and we obtain: 


dw v dv ln(P 2 /P, ) 

dz r dr z, - z. 





This equation is a first order differential equation which is solved by numerical computer 


techniques by the method of finite differencing. To do this, the equation is expressed in 
finite difference form as follows: 


w 


i+i.j 


- w 




Az 


v i,j v i,j+i " v i,j _ ln(P i+ ,/Pi) 

— — Wj j 

fj A f z i+ 1 “ z i 


(29) 


where Az and Ar are chosen grid spacings for altitude and latitude, respectively, and i,j are 
level indices. The above equation is readily solved for w i+ , : if w ; , is known. Hence the 
boundary condition was specified to be Wj j = 0 at an altitude of 15 km to approximate 
the real physical situation. The values of the v’s, P’s and z's were known from previous 
calculation of the meridional speed and from the input pressure contour data. 

t. INPUT DATA AND NUMERICAL COMPUTATION OF WIND SPEEDS 

The algorithm for calculating the zonal, meridional and vertical wind speeds according 
to the polar vortex model was coded in FORTRAN and run on the IBM S/95 Computing 
System at Goddard Space Flight Center (GSFC). The input data for the model were ex- 
tracted from pressure contour maps that exhibited nearly circular vortices around the 
North Pole, examples of which are given in Figures 4 and 5 for the 2 mb and 0.4 mb sur- 
faces, respectively. Each chart represents a weekly or monthly average of the measured 
data. Figure 6 shows a summary of the input data obtained from the charts for the 0.4, 

2, 5, 10, 30, 50 and 100 millibar surfaces, measured during October 1967. Each curve 
represents the variation of the height with latitude of each isobaric surface. It was found 
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Figure 4. Contour Map of the 2— Millibar Pressure Surface, Showing Weekly Averages 
for Oct. 11, 1967. (ESSA Technical Report WB 11.) 
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Figure 5. Contour Map of the 0.4-Millibar Pressure Surface. Showing Weekly Averages 
for Oct. 1 1. 1967. (ESS A Technical Report WB 12.) 




Figure 6. The Heights of Various Isobaric Surfaces as Functions of 
Latitude, Extracted from Contour Maps for Oct. 1967 


17 


that the late fall/winter months was the time of the year when there was the best likeli- 
hood for finding well developed circular vorticies. 

As described in Section C above, data values of Rj n , R r ef,n ant * z ref,n were ‘Ob- 
tained from the n 1 ^ isobaric contour. R re f jn and z re f n were chosen to be the end point 
of the n*b curve. These values were used to compute the best values of the vortex 
strength K n and vortex radius A n at the average height of that contour. A fifth degree 
polynomial was then fitted by least squares to the K’s for the several isobars. A similar 
fit was made for the A’s. Thus K and A were obtained as functions of altitude z, the 
plots of which are shown in Figures 7 and 8. Equation (15) was then used to calculate 
the zonal wind speed u at grid points (R, z) in the latitude range 90°-45°N and the alti- 
tude range 15-60km. A contour plot of the resulting zonal winds is shown in Figure 9. 

The kinematic Prandtl eddy viscosity was then calculated at each point according to 
Equation (18), using a mixing length of 10 meters, which is close to estimates made by 
Holton (1972). The resulting kinematic eddy viscosity contour distribution is shown in 
Figure 10. Figure 11 shows a cut of the Prandtl viscosity profile at 40°N latitude, and 
also values used by other investigators. The profile at 80°N is shown in Figure 12. It 
is seen that the profile derived from this model at a given latitude exliibits a wavelike 
variation with altitude. The magnitudes of the eddy viscosities are of the order of 10 
ineter 2 /sec in agreement with values used by Cunnold, et ai., (1975) and by Justus (1973). 
Tiie meridional wind speed v at each grid point is then determined from Equation (19), 
resulting in the contour plot shown in Figure 13. 





ZONAL WIND IN METER5/5EC0ND 



Figure 9. Contour Distribution of Zonal Wind Derived from Oct. 1967 Data 


V 




21 



m o c -h 


A 

L 

T 

I 


1 

N 

K 

M 


60 r 
55 


50h 


45 ^- 


i r 


I 


4Qi 


35 


30' 


25 


20 


15h* 


10 L - + 1 

0 2 4 




LEGEND 

j POLAR VORTEX MODEL 

* FROM CUNNOLD ET AL. 1 
A FROM JUSTUS 


i 


■'K, 


L 


i 


— 

X 




7 ^ 










1 




j 






6 8 10 12 M 

EDDY VISCOSITY IN M**2/5EC 


18 20 


Figure 1 1. Prandtl Kinematic Eddy Viscosity Profile at 40°N Latitude from 
the Polar Vortex Wind Model Using Data from Oct. 1967. Also shown are 
values of eddy viscosity used by other investigators. 
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EDDY VISCOSITY IN M**2/5EC 


Figure 12. Prandtl Kinematic F.ddy Viscosity Profile at 80°N Latitude from the 
Polar Vortex Wind Model Using Oct. 19o7 Data 
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The vertical winds w at the prill points were calculated by the recursion relation 
Equation (24), imposing the condition that the vertical wind component is zero at the 
lower boundary at 15 km. This condition is a reasonable assumption since it is known tnat 
vertical motion is negligible at that height. The resulting vertical wind contours are given 
in Figure 14. 

Ci. DISCUSSION OF RISULTS AND CONCLUSIONS 

Using measured pressure contour data as input, the polar vortex wind model is used 
to calculate the Prandtl eddy viscosity distribution and the associated stratospheric wind 
components in the latitude range 40°N-45°N and at altitudes between 15 and 60km. The 
magnitudes of the kinematic eddy viscosity, derived by assuming an experimentally reason- 
able mixing length, were found to approximately agree with values used by Wofsy and 
McElroy (1473) in duplicating the observed distribution of methane in a one-dimensional 
chemical-diffusive model, and by Cunnold, et al., (1475) and Justus (1973). However, 
unlike the others, the eddy viscosity distribution obtained and used by the present model 
exhibits pronounced latitude-altitude dependence with a wavelike vertical structure in any 
given latitude plane. 

The magnitudes of the zonal, meridional and vertical winds derived by this model were 
of the order of meters, centimeters and millimeters, respectively, in agreement with results 
published by the other investigators (Cunnold, et al.. and Vincent, 1463). It was found 
that for cyclonic zonal llow in the Northern Hemisphere, the flow direction in a meridional 
plane is radially toward the pole al lower altitudes, increasing upward at intermediate 
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Figure 14. Vertical Wind Contour l 


levels, and radially outward at higher altitudes, as expected from fluid-dynamical 
considerations. 

The present model is capable of utilizing experimental satellite data for deriving the 
corresponding eddy viscosity and wind distributions, in a semi-empirical manner. Further- 
more, by being a purely fluid-dynamical < ?del, it avoids the problem of circular reason- 
ing characteristic of the Wofsy model. The senior author has always been ill at ease with 
the reasoning which derives eddy diffusivity by matching constituent profiles in photo- 
chemical mod e Is, then using the eddy diffusivity thus derived for transport in other chem- 
ical kinetic models. Also, the empirical nature of our model permits a study of variations 
of eddy transport with latitude and time of year, and from one year to another. 
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